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Abstract 

We investigate the time-evolution and steady states of the stochastic susceptible-infected- 
recovered-susceptible (SIRS) epidemic model on one- and two- dimensional lattices. We compare 
the behavior of this system, obtained from computer simulations, with those obtained from the 
mean- field approximation (MFA) and pair-approximation (PA). The former (latter) approximates 
higher order moments in terms of first (second) order ones. We find that the PA gives consistently 
better results than the MFA. In one dimension the improvement is even qualitative. 

PACS numbers: 87.23. Ge,05.70.Ln 
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I. INTRODUCTION 



The mathematical modehng of the spread of epidemics is a subject of continuing theo- 
retical and practical interest |^ . This is enhanced by the fact that the same or similar 
models are used for describing other phenomena such as plant and animal dispersal, and 
successional dynamics in ecology 

The level of description provided by a model can be purely macroscopic and deterministic 
or individual and stochastic j^. In the first case one uses (partial-) differential equations 
to describe the time evolution of different subpopulations; e.g., susceptible, infectious and 
recovered. In the second case one typically uses stochastic dynamics on a lattice (or more 
general graphs) where the variables at each node represent the state of an individual or a 
small spatial region. The time evolution of these variables is stochastic, e.g., an infected 
individual at site i has a certain probability per unit time (rate) A to infect a susceptible 
individual at a neighboring site j. These systems fall into the category of what mathemati- 
cians call interacting particle systems 0, Q| and physicists call stochastic lattice gases - 
systems of great interest also in the study of equilibrium phase transitions, phase segregation 
kinetics, etc., fields very different from epidemiology and ecology. 

The connection between these modes of description and various intermediate ones has 
been i.ves.ig.ed extensively ,n .eeent yea., e.g.. see flafly QQ. M*e.a.ca..y 
this involves the use of the so called hydro-dynamical scaling limit. This uses a rigorous 
separation of space and time scales to derive deterministic macroscopic equations from the 
microscopic dynamics of stochastic lattice systems. Other approaches are based on more 



heuristic methods such as the mean field approximation (MFA) and improvement thereof 
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The present work falls in the latter category. We apply a pair approximation (PA) scheme 
to a microscopic stochastic epidemic model in which individuals recovered from an infection 
enjoy a period of immunity before again becoming susceptible at a rate 7: the SIRS model. 
The PA approximation was used by Durrett and Levin for the simpler susceptible- 
infected-susceptible (SIS) model where recovered individuals immediately become suscep- 
tible again. They compared the results of the PA and MFA with those of the stochastic 
SIS model and found that the PA gave a quantitative improvement over the MFA. Here 
we consider the general SIRS model. We obtain the behavior of the stochastic model from 
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extensive computer simulations. We then solve the PA and MFA models analytically for the 
stationary state and numerically for the time dependent case. We find that the PA gives con- 
siderably better agreement with the simulations than the MFA both for the time evolution 
and for the steady state. For the latter the PA reproduces the qualitative difference between 
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the one and higher dimensional phase diagram of this model found in Ref. 
This is reminiscent of the relation between the MFA and the Bethe-Peierls approximation 
(which the PA closely resembles) for equilibrium lattice systems 



II. THE STOCHASTIC SIRS MODEL 

We first recall the stochastic lattice model of the SIRS epidemic process j31|. A site x 
of a (i- dimensional lattice can be occupied by an individual in a state of S (healthy and 
susceptible), I (infected), or R (recovered, i.e., healthy and immune). The system evolves 
according to the following transition rates, 

S I at rate An(x), (1) 
I R at rate 6, 
R S at rate 7, 

where n{x) is the number of infected (nearest) neighbors of x, A is the infection rate, S is the 
recovery rate and 7 is the rate at which immunization ceases. The limit 7 ^ 00 corresponds 
to the case where a recovered site passes instantaneously through the state R; this is the 
SIS model, also known as the contact process. We shall choose time units in which 6=1. 

One can obtain some rigorous qualitative information about this and related models via 
p..obabn,stic app..o..hes such as tKose used in .nt.act.g pa.t,c>e systems 00^0. 



Of particular interest is the behavior of the stationary state on an infinite lattice which is a 
good approximation for the quasi-steady state behavior of large systems: see Appendix |El 
. This information is encoded in the phase diagram of the stationary state which depends 
on the infection rate A, the recovery rate 7 and the topology of the lattice. For small A, 
the only stationary state is one in which all sites are in the susceptible (disease-free) state 
while for large A there is (for the infinite system) also a stationary state containing non zero 
fraction of I and R individuals. 

The critical infection rate Ac(7) is defined as the smallest value of A, for a given 7, above 



which the infection can persist forever. For the SIS or contact process (7 = 00), the critical 
infection value is known with high accuracy, Ac(oo) ~ 1.6489 in d = 1 and Ac(oo) ^ 0.4122 
in = 2 B 3. Considerably less is known about the phase dia Efram of the SIRS model. 
Interestingly there is a qualitative difference in the behavior of Xdl) in one and in higher 
dimension when 7-^0. It has been shown that lim^^QXc{'~f) = Ac(0) is finite when d > 2 
while Ac(0) = 00 when d = 
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To go beyond qualitative results we need to carry out simulation or make some approxi- 
mations. This is the subject of the rest of the paper. 

III. THE PAIR APPROXIMATION 



The time evolution of the single site probabilities in the stochastic SIRS epidemic process 
can be written in the following form. 



dt 

dt 

dPt{R.) 
dt 



-A Pt{S.Jy)+lPt{R.). 

A Pt{S.,Iy)-Pt{Q, 

Pt{h)-iPt{R.), 



(2a) 
(2b) 
(2c) 



Here M{x) is the neighborhood (nearest neighbor sites) of a site x, Pt{oix) is the probability 
of having a state a at site x at time t and Pt{o.xi (3y) is the joint probability to have state a 
at site X and state [3 at site at time t. We always have Pt{Sx) + Pt{Ix) + Pt{Rx) = 1- 

Eqs. fl2ap - (Pcj) are, as is usual for moment equations, not a closed system. One can extend 
them by including equations for the time evolution of Pt{Sx, ly) which in turn involve higher 
moments of the spatial correlations. This leads to an infinite hierarchy. To solve such 
a hierarchy one usually resorts to some approximation scheme which expresses the higher 
order moments in terms of the lower order ones and truncates the equations at some point ; 
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this is referred to as the moment closure method ^ , 
Both the MFA and PA are such schemes. In the MFA Eqs. (|2a| ) - (Pcj) are closed by assuming 
that Pt{Sx,Iy) = Pt{Sx)Pt{Iy), i.e., it neglects correlations between different sites. This 
leads to a pair of coupled equations which have been studied in jsil. In the PA scheme 
Pt{<yx) and Pt{ax,Py) are kept as unknowns while the higher-order moments are expressed, 
via some appropriate approximation, in terms of these quantities. 



To carry out the PA we complement Eq. Q by equations for the second moments 
Pt{ax, Py) for nearest neighbor sites x and y based on the transition rule that we have 
described in Eq. ((T)). These are 



dPt{Sx, ly 
Jt 



-- ^Pt{RxJy)-{\ + l)Pt{SxJy)+ Yl ^Pt{Sx,Sy,I^) (3a) 

~ ^Pti^w, Sx, ly) , 

"^^'^^^l^y^ = Pt(Sx, ly) + iPt{Rx. Ry) - iPt{Sx, Ry) - Yl ^y^^ (3b) 

"^^'^^r^'^ = -{l + l)Pt{Rx,Iy)+Pt{Ix,Iy)+ Y ^Pt{Rx.Sy,I^) (3C) 

where J\f^{y) is the set of nearest neighbor sites of y excluding the site x. Ptidx, Py,Xw) is 
the joint probability to have state a at site x, state j3 at site y and state x site w at time 
t. For a derivation of Eq. (jS)) see Appendix 1X1 

To close the system Q and Q and derive a set of autonomous equations for Pt{ax) 
and Pt{ax, Py) we approximate the triad joint probability Pt{ax, Py, Xw) for x and w nearest 
neighbors of y, by the product of two pair probabilities Pt{ax, Py) and Pt{Py, Xw) divided by 
the probability Pt{Py) 14, 3, 12, 3, i.e., we set 



p. a Pt{ax,Py)Pt{Py,Xw) . . 

r't[ax,Py,Xw) B77r\ 

^t[Py) 

Note that we have made use here of the structure of the hypercubic lattice. In such lattices, 
three adjacent sites, x,y,w can not form a triangle but form only linear chains. This is 
not so in other lattices, e.g., the triangular lattice, where other configurations need also be 
considered. 

While there are other choices for a PA, the approximation in Eq. (jlj) allows one to 
get the steady state solutions analytically. With other pair approximations jssi, one has to 
solve the resulting differential equations numerically, making it impossible to obtain analytic 
expressions for the critical curve. 

To actually carry out computations with the PA we will assume from now on that our sys- 
tem is spatially uniform. The site x, in Eqs. and (jH]), can now be chosen to be the origin. 
We also define: Pt{S,I) = \ T.y&M{x)Pt{Sx,Iy), Pt{a,p,x) = E«,eAr-(y) /^y 
where z = 2d is the number of nearest neighbors of a site in the d-dimensional cubic lat- 
tice. The truncated equations for the PA-SIRS can now be written, by using the exact 



Eqs. (121)- Q and the approximate Eq. (jH), as a closed set of five coupled equations, 
dPtil) 



dt 
dPtjR) 
dt 



zXPt{S,I)-Pt{I) (5a) 
Pt{I)-^Pt{R) (5b) 



dPt{S,R) ^(r,lTy\ nfjD T\ OP/Q jyw - ^)^PtiS,I)Pt{S,R) 

= P,{S, I) + 7(P.(i?) - Pt{R, I) - 2P,iS, R)) - i_p^(^)_p^(j) (5c) 

^ ^ -(2 . ,WR, /) + P.(/) - P.iS, /) + ^y^^t^'!^;-^^^^^^ (5d) 

^^^^ = 7Pt(P, /) - (A + l)Pt(5, /) (5e) 

Note that we always have Pi(a) = Pt(a, S) + Pt(a, /) + Pt(«, P) which determines Pt{I, I) 
and Pt(5,S). 

In the limit 7 ^ 00, Pt{R) and Pt{R, a) as well as their time derivatives will go to zero. 
This yields 7Pt(P) = Pt{I) and jPt{R,I) = Pt{I) - Pt{S,I) m^. In this limit Eq. © 
reduces to the PA equations of the SIS considered in jl6| . 

"^^'^^^ = zXPt{S,I)-Pt{I), (6a) 



dt 



^^1^ = P,(/) - (A + miS, I) + i^jp^ff '^ (1 - W - 2^*(^' ^))- (6b) 

As already noted the MFA approximates the joint probability Pt{S, I) in Eq. ()5a|l by the 
product, PfjS, I) = Pt{S)Pt{I). This leads to the closed set of MFA of equations for the 

SIRS lail, 

= -zXPt{S)Pt{I)+^Pt{R) (7a) 

dPt{I) 



dt 
dPt{R) 



zXPt{S)Pt{I) - Pt{I) (7b) 
Pt{I)~iPt{R) (7c) 



dt 

For 7 ^ 00, 7Pt(P) Pt{I) and Pt{S) 1 - Pt{I). Eq. © then reduces to the MFA for 
the SIS. 



IV. STATIONARY SOLUTIONS OF THE PA-SIRS MODEL 

Let us first consider the steady state solutions of the PA-SIS obtained by setting the l.h.s 

^ n 

of Eq. (jo)) equal to zero |1^] . This gives for the critical value of the PA-SIS epidemic process 
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Ac(oo) = 1/(2; — 1). For A < Ac(oo), both Pt{I) and Pt{S,I) as t — > 00 for all initial 
states. When A > Xc{oo) there is, in addition to the disease- free state corresponding to 
P(/) = 0, also a stationary state consisting of a finite fraction of infected individuals: 

Pis, I) = P{I)/{zX) (8a) 

It is these non-zero steady states which are approached as t — >■ 00 when starting from any 
initial state with Pq{I) > 0. 

The steady state solutions of the PA-SIRS system is obtained by setting the l.h.s. of 
Eq. (0) equal to zero. Setting x = P{I) this yields, 

P{R) = x/7 (9a) 



P{S,I) = x/{zX) (9b) 

P{S,R) = (9c) 

= ^-^^'f^' (9d) 

1 I 1 



where P{a,f3) are the approximate probabilities for having states a and f3 on neighboring 
sites. After further simplifications, we find that x has to satisfy the cubic equation, 

x{aix'^ + a2X + as) = (10) 

Both the derivation of Eq. (jlUj) and the explicit expressions for ai, a2 and as functions of 
A and 7 are given in Appendix FBI 

The root x = corresponds to the all healthy steady state, which is always a solution. 
The critical curve Ac (7) is determined by the existence of a root of Eq. (fTUI) such that x 
and all other stationary probabilities, are strictly positive. It turns out that this strictly 
positive root is unique. Thus when A < Ac(7), a; = is the only steady state solution. For 
A > Ac(7), there is also a steady state in which the infection is endemic: P{I) = 'yP{R) = x 
and P{S) = 1 — (1 + l/7)x, see Appendix |Bl 

The critical curve Ac(7) is obtained in Appendix |B| It is given by the equation. 



As 7 ^ oo, Ac(oo) = (2d — the critical point of the PA-SIS epidemic process. On 

the other hand, as 7 approaches zero, the critical curve shows different behavior depending 
on the dimension of the lattice: Ac(0) diverges to infinity for d = 1, while Ac(0) is finite 
for d > 2. The PA thus reproduces the qualitative difference between the one and higher 
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dimensional phase diagram of the SIRS model found in Ref. 

The MFA, Eq. ((Tj), yields the mean field critical value, = 1/z independent of 7. In 
the coexistence region A > X^^^ the mean field stationary states are P{I) = ^P{R) = 
and P{S) = ^. 

Both the steady state and critical value of the MFA and PA fail to correctly represent 
the results of the stochastic SIRS process for small 7: see Figs. nand|21 Note in particular 
that P{S) of the stochastic SIRS process is considerably larger than that of the MFA or PA 
for large A and small 7. This is due to the fact that the susceptible sites can be surrounded 
by recovered ones and thus protected from contacting infected ones in the stochastic case. 



V. COMPARISON OF THE STOCHASTIC, THE PA AND MFA STEADY 
STATES. 

We compare in Figs. El- El the steady state values of -P(a) and P{a,P) obtained from the 
MFA and PA with the results from the stochastic SIRS process as a function of A at fixed 
values of 7. Clearly the PA gives results closer to those obtained from the stochastic model. 
For the methods used to obtain the steady state results from the numerical simulation, see 
appendix. lEl 

Figs. El and show that both the MFA and PA overestimate P{I) as well as P{a,I), 
a = S,R. This is due to the strong tendency of infected sites in the stochastic model to 
cluster into localized islands, reducing the contacts between S and /. This is partially taken 
into account by the PA as seen by the behavior of P{S, I) and P{I, I) in Figs. El and This 
clustering effect is also observed in the stochastic SIS process [igI . It is more pronounced in 
one dimension. 

Note that P{S, I) becomes zero both at A < Ac(7), when P{I) = 0, and at A = 00 when 
P{S) = 0, reaching a peak at a positive value of A which depends on 7. For large values of 
7, the steady state values of -P(a) and P{a, (3) obtained from the PA, or the MFA agree well 
with the numerical simulation, away from the critical Ac(7). Moreover the PA yields steady 
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state curves remarkably similar to those from the numerical simulation, see Figs. E] and IHl 



VI. LINEAR STABILITY ANALYSIS OF THE PAIR APPROXIMATION 

To study the stability of the stationary PA state, Eq. (jSJ is linearized about the steady 
state values jsil, see Appendix O This leads to the study of the roots of the characteristic 
fifth order polynomial P^i^), obtained from |y4 — ^/| = where A is the Jacobian of the 
linearized PA-SIRS system. If ReC, < 0, the solution of the linearized equation is stable, i.e., 
a small perturbation around the steady state will decay back to the steady state. We used 
the Routh-Hurwitz conditions ,31] to obtain the sign of the real part of eigenvalues of the 
Jacobian. As expected, the positive steady state solution is stable for A > Xdl)- The zero 
steady state solution is stable for A < Ac (7) and unstable for A > Ac (7). 

The eigenvalues of P^i^) have non-zero imaginary parts in some regions of the parameter 
space. In such regions the PA-SIRS system in Eq. ^ will converge to the steady state in a 
damped oscillatory manner. Such oscillations are seen in Fig. [3 and |H1 



VII. TIME DEPENDENT BEHAVIOR 



To study the time evolution of an epidemic following an initial infection of a healthy 
population we performed dynamical Monte Carlo simulations js^ as well as solutions of 
Eqs. (0) and ((Tj). For the stochastic evolution we started with infected sites placed either 
randomly or in a cluster and followed the time evolution averaged over 10^ realizations of 
the SIRS process. To obtain the time evolution of the MFA and PA we solved Eq. ((7j) and 
Eq. (jS)) numerically by using a 4th-order Runge-Kutta method. We plot the results in Figs. [7| 
andEl 

To set the unit of time of the simulation we started with a fully infected state, Po{I) = 1 
and A = and obtained the exponentially decaying pattern of Pt{I)- We then set the slope 
(death rate) of the graph, logPt{I) vs t, from the numerical simulation equal to those from 
the MFA and PA. 

Starting with a small value of Po{I), Pt{I) displays an initial "exponential" growth in 
both the MFA and PA. Similar growth patterns are observed in all Pt{a,I), a = S,I,R. 
This is explained by the initially abundantly available susceptible population. Once the 
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susceptible population is reduced, the infected population reaches a maximum and then 
decreases to the steady state endemic level. Note the damped oscillatory pattern in Figs. [7| 
andlHlfor this choice of the parameters (A,7). 

The numerical simulation of the stochastic time evolution does not show the pronounced 
growth patterns of the PA and MFA when the initial fraction of infected sites is small, as seen 
in Fig.|H| The formation of clusters of infected sites makes the infected population grow more 
slowly in the stochastic model. When the initial fraction of infected population increases to 
more than one percent the stochastic model shows significant change in its growth pattern, 
becoming similar to the PA and MFA. If however the same fraction of infected sites are 
initially placed in a single cluster the stochastic epidemic process exhibits slower growth 
patterns, similar to those starting with a small fraction of initially infected sites. These 
studies confirm that the clustering of infected sites in the stochastic model reduces both 
the speed of growth and the maximum fraction of infected sites. In realistic situations the 
population is not well mixed so we would expect growth patterns more similar to that of 
the stochastic epidemic model, starting with a fraction of infected sites initially placed in a 
single cluster. 

VIII. SUMMARY 

We investigated the stochastic SIRS epidemic process and compared the results with those 
obtained from the deterministic MFA and PA. These approximations close the hierarchy of 
dynamical equations by expressing the higher order moments in terms of the lower order 
ones. The PA is found to improve over the MFA both for the stationary and for the time 
dependent states. The time evolution of the system shows damped oscillatory behavior in 
some parameter ranges. 
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APPENDIX A: DERIVATION OF DIFFERENTIAL EQUATION FOR Pt{S^Jy) 

Eq. ()3a|) is derived by considering all transitions leaving or entering the pair configuration 
{Sx,Iy). We list them as follows: A pair {Rx.Iy) changes to a pair {Sx.Iy) with a rate 7. 
A pair {S^, ly) changes to a pair (J^., ly) with a rate A and also changes to a pair {S^, Ry) 
with a rate 1. A triad configuration {Sx-, Sy, 1^) transits to a triad {Sx, ly, Iw) with a rate A 
such that a pair configuration (5*^., Sy) is changed to (5^., ly). A triad (/^, S^;, J^) changes to 
a triad (J^, Ix, ly) with a rate A. The equations for Pt{Sx, Ry) and Pt{Rx, ly) in Eq. Q can 
be obtained in a similar way. The relation Pt{ax) = Pt{cix, C(y) + Pt{c(x, Py) + -Pt(ttx5 Xy) can 
be used to obtain the other joint probabilities Pt{ax,Py) which are not shown in Eq. (jS)). 

APPENDIX B: DERIVATION OF EQ. ^ 

The steady states in Eq. © are obtained by setting the l.h.s. of Eq. (jHaj) - (|5djl equal to 
zero. In addition we set Eq. (j3e|) equal to zero and replace a single site and joint probabil- 
ities with the steady states in Eq. Q. After simplifications, we obtain Eq. ()10|) with the 
coefficients, 

ai = -f^{z^{z - 1)A -z}+ -f^{z{2z^ - 2z - 1)A - 2z - 1} (Bl) 

+ -f{2z{z^ -z-l)\-2z~l} + z{{z^ - 2 - 1)A - 1} 
02 = z-^^-i^{z + 1 - 2z{z - 1)A} + 7{z + 3 - (3^2 - 4z - 1)A} + ^ + 1 - {2z^ - 3^ - 1)a| 
as = A'{7{-1 + - 1)} - 1 + \{z - 2)}. 

The critical curve Ac (7) is given by setting 03 = 0. Only for A > Ac (7) does the quadratic 
factor of Eq. have a positive root. 
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APPENDIX C: THE JACOBIAN OF THE LINEARIZED PA-SIRS 



The Jacobian of the hnearized PA-SIRS is written, 



A 



-7 


'K2K0 



1 

-1 
'K2K0 





zX 



P(IS) 



where Kq = 1 + 2 



7-7^2 -^2 1 - 

P(IS) ^ ^ (z-l)XP{IS)P{SR) ^ _ 








P{SR) 



-27- 



P{SR) 



P{SR) 






7 

-7 

-7 - 



2/ 



P{SR) 



1~P(R)-P(I) 



{z-l)XP{IS)P{SR) 
(l_P(ij)_P(/))2 : 



and K3 = {z- 2)A - 



1 



+ 



P{IS) ' P{SR) ' 



APPENDIX D: LINEAR STABILITY ANALYSIS OF THE MF-SIRS 

The Jacobian matrix B of hnearized MF-SIRS is given by [sij] 

^^l ~XzP{I)-^ _AzP(S)-7 
\^ XzP{I) XzP{S) - 1 

The characteristic polynomial of the second order, ^2(0 = 'C^ + o^ii + 02 = 0, is obtained 
from \B - n\ = 0. 

rn 

The necessary and sufficient (Routh-Hurwitz) conditions 31] for ReC, < is 02 > and 
ai > 0. In the coexistence region where zX > 1, 02 = 7(zA — 1) > and ai = ;^(7 + 2;A) > 
for all 7 > 0. In the no-coexistence region where zX < 1, 02 = 7(1 — zX) > and 
ai = 7 + (1 — zX) > for all 7 > 0. Both in the coexistence and no-coexistence region, the 
real part of the eigenvalues is negative and thus the mean field steady states are stable. 

Now we turn our attention to the oscillatory behavior. The eigenvalues of the character- 
istic polynomial ^2(0 is given by, 

-72 - ;z7A ± v/(27-22A2 - 2zX{'y^ + 2z-fX + 2) + 7^ + 47^ + 87 + 4)7 

= WiTT) ^^'^ 

In the range of A_(7) < A (7) < A+(7) the imaginary part of the eigenvalues is non-zero: 
■^±(7) = ^^"^'^"'"'^ ^7^^^"'"'^'' ^ • ^^^^ range of A, the steady states correspond to the stable 
spiral and the system converges to the steady state in damped oscillatory pattern. Even 
in the damped oscillatory region, any oscillation is hardly visible in the large 7 limit and 
becomes noticeable only in small 7 limit. 
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APPENDIX E: MONTE CARLO SIMULATION 



The numerical simulations described here used lattices with periodic boundary conditions. 
In one dimension, rings of 5000 < N < 15000 sites were used. In two dimensions, torii of 
50^ < N < 200^ sites were employed. 

To obtain the steady state of the SIRS process a random initial configuration of susceptible 
and infected sites is evolved according to the transition rates in Eq. (^. In practice a site 
is randomly chosen and a random number (g [0, 1]) is also chosen: if it is greater than the 
given transition probability for that site, which is equal to the ratex At, its state is updated: 
At is chosen to be so small that transition probability is not greater than 1 for a range of 
i^^l) ^Is^- Otherwise its state remains the same. 

For a finite system the only true stationary state of the SIRS process is the absorbing 
state corresponding to P{S) = 1, P{I) = P{R) = 0. To learn about the active state from 
simulations of a finite system we study the quasi- stationary state. These are determined 
from averages over the surviving representatives of 10^-10^ independent realizations of the 
SIRS process with the same parameter (A, 7), beginning with random initial distribution of 
the /'s. Surviving sample averages converge to stationary values as — 00. To obtain the 
steady states and critical curve we extrapolated quasi-stationary values of finite systems to 
those of the infinite system. 

The finite size scaling theory ^] can be used to obtain the critical curve X^il)- We can 
assume a scaling function of the surviving probability: Pt(/) ~ r^/'^ii/((A - AJt^/'^")- At 
criticality, A = Xdl), the survival probability of the infection, starting from a single infected 
site, has a power law behavior in time. In the subcritical region, it decays exponentially 
while in the supercritical region it reaches non-zero steady state in a short time. The power 
law behavior of the survival probability at criticality enables one to extract the critical 
curve X^il) fro^i the time evolution data of the SIRS process. This dynamical Monte Carlo 
simulation is reliable when the system size is sufficiently large so that the evolution of the 
system is approximately confined, for the duration of the simulation to a region smaller 
than the size of the system Q]. However we found that this surviving probability oscillates 
wildly when 7 is small. Because of this the dynamical Monte Carlo method is not used to 
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determine the critical curve near 7 = 0. 
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FIG. 1: Phase diagram of the SIRS process in two dimensions. The coexistence phase of S-I-R and 
the no-coexistence phase are separated by the critical curve from the simulation (open circles with 
dotted line for eye-guidance) , the PA (thick solid line) and the MFA (long dashed line) . The critical 
curve is obtained on periodic square lattice of different sizes N from simulations extrapolated to 
infinite system : iV = 50^, 70^, 100^, 150^, 200^. 
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FIG. 2: Phase diagram of the SIRS process in one dimension. The critical curve from nu- 
merical simulations of ring lattice of different sizes N is extrapolated to infinite system: N = 
5000, 7000, 10000, 15000. The same symbols are used as in the Fig. [H 
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FIG. 3: First and second order moments of the steady state SIRS in two dimensions at 7 = 0.2. 
The steady-state values of the density of infection in Fig.|3fa) and the second moments in Fig.lHfb)- 
(f) are drawn from the numerical simulation (open circle with dotted line for eye-guidance), the 
PA (thick solid line), and the MFA (long-dashed line). For the numerical simulation we used a 
system of size = 100^. 
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FIG. 4: First and second order moments of the steady state SIRS in two dimension at 7 = 2. The 
same symbols are used as in the Fig. |3I 
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FIG. 5: First and second order moments of the steady state SIRS process in one dimension at 
7 = 1. The same symbols are used as in the Fig. |3I 
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FIG. 7: Time-evolution of the first and the second order moments of the SIRS process in two 
dimensions. All sub-graphs are from numerical simulations (open circles), the PA (solid line), 
and the MFA (dashed line) at 7 = 0.2 and A = 2. A periodic square lattice of N = 10^ sites is 
used in the numerical simulations averaged over 10^ — 10^ realizations starting with random initial 
distribution with 1 percent of infected sites. 
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FIG. 8: Time-evolution of fraction of infected sites of the SIRS process in two dimensions at 7 = 0.2 
and A = 2. A periodic square lattice of TV = 100^ is used in numerical simulation averaged over 
10^ — 10^ realizations. Main: Simulation starts with 1% of infected sites placed either randomly 
(filled circles) or in a single cluster (open circles) on a lattice. Both the PA and MFA takes an initial 
value 0.01 for Po{I)- Inset: Simulation starts with different fractions of infected sites randomly 
placed in a lattice: 0.1, 1 and 5 percents of the system. 
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